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ABS TRACT 


A computer simulation model is presented of soil water movement and 
uptake by plant roots including the computer software package which 
implements the model. The flow of soil water is assumed to satisfy 
Darcy's law and is modeled with a nonlinear parabolic partial differential 
equation. Soil-water uptake by the root system is accounted for in the 
model by a volumetric sink term. The numerical implementation of the 
model uses one dimensional linear elements and an implicit time step. 

The model is designed to be sufficiently versatile that it can be applied 
to a wide range of soils and function with available daiiy meteorological 
data. 


Efficient and reliable computer software is a prerequisite to 
implementing an understanding of soil physical processes into operational 
prediction schemes. Another requirement is more routine methods of 
specifying soil water characteristics. The bulletin is designed to 
provide the technical documentation of a computer software package that 
was developed with a view to contributing in these two areas. It 
provides the necessary technical material so that the mathematical 
formulation, the numerical techniques and the computer implementation 
can be assessed. 


RESUME 


Les auteurs présententiiun modéle de simulation mécanographique du 
transport de l'eau dans le sol et de son absorption par les racines 
des plantes, ainsi que l'ensemble de programmes qui appliquent le 
modéle. La migration de l'eau est censée satisfaire la loi de Darcy 
et est modélisée par une équation différentielle partielle, parabolique 
et non linéaire. L'absorption de 1'eau du sod par le systéme racinaire 
est prise en compte dans le modéle par une variable de puits volumétrique. 
L'application numérique du modéle fait appel a des éléments linéaires 
unidimensionels et a une technique implicite par avancement de temps. 
Le modéle est congu de facgon suffisamment souple pour pouvoir s'appliquer 
a une grande diversité de sols et fonctionner 4 méme les données météo- 
rologiques quotidiennes disponibles. 


L'utilisation d'un logiciel efficace et fiable est la condition 
préalable a la compréhension des processus physiques du sol et a4 leur 
application dans des plans prévisionnels d'exploitation. Une autre 
condition consiste a utiliser des méthodes plus courantes de détermi- 
nation des caractéristiques de l'eau du sol. Le bulletin a pour objet 
de fournir la fiche technique d'un ensemble de programmes élaborés en 
vue d'enrichir nos connaissances dans ces deux domaines. I1 fournit 
le matériel technique nécessaire pour évaluer la formulation mathéma- 
tique, les techniques numériques et la mise en forme mécanographique. 


INTRODUCTION 


Mathematical models can provide a useful synthesis of soil science, 
plant science and agrometeorology (Philip, 1969; Hillel, 1971). 
Empirical statistical approaches are often best for prediction for the 
average situation, but models which simulate processes can provide 
better insight into one year out of ten occurrences where there would 
not be adequate data for a valid statistical prediction. In addition, 
the largest benefit from mathematical modeling generally results from 
the analysis of abnormal situations. It is felt that this fact provides 
further justification for incorporating relevant physical and biological 
theory in models. 


In order to solve many of the existing soil water simulation 
models, detailed soil, crop and climatic data are required (see eg. 
Nimah and Hanks, 1973, Feddes et al., 1974; van Bavel and Ahmed, 1976; 
van der Ploeg et al., 1978). Consequently the models have generally 
been site specific and this has limited their widespread use. 


The objective of the present work is to describe a computer 
simulation model of soil water movement and water uptake by plant 
roots. The model is designed to be sufficiently versatile that it can 
be applied to a wide range of soils and function with available meteorolo- 
gical data. Simulations and predictions of the soil water regime in 
the rooting zone and actual transpiration for a growing crop can be 
carried out. Since the ratio of actual to potential transpiration can 
be related to the growth rate of a crop (Hanks, 1974), the model can be 
considered as a tool for the study of soil water effects on yield. 


The flow of soil water is assumed to satisfy Darcy's law and is 
modeled with a one-dimensional non-linear parabolic partial differential 
equation. Soil water uptake by the root system is accounted for in the 
model by a volumetric sink term (Feddes et al., 1974). The surface 
boundary condition is determined from daily precipitation, potential 
evapotranspiration and a specification of crop cover. The basal boundary 
condition can be specified either as a flux or a soil water concentration. 
The option of specifying a water table is also incorporated. 


The computer program is written with the objective of providing 
efficient and stable numerical algorithms that make it feasible to 
carry out numerical simulations for the entire growing season. Crop 
information relating to ground cover and root distribution can be 
specified as a function of crop development stages. Standardized 
functions are used for soil water relationships (Clapp and Hornberger, 
1978), and up to four distinct soil layers can be simulated. Only 
minimal changes are required to simulate different soil types. 


In the following, the FORTRAN computer program is described, 
including details required to set it up and run it. The main computer 
program and subroutines are discussed in sufficient detail so that 
it should be clear where changes are required to adapt it to a range 
of possible simulation problems. It is assumed that the implementation 
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of a diffusion-based soil water simulation model would require a 
knowledge of soil physics sufficient for understanding the implications 
of specifying soil physical properties and adequate knowledge of 
numerical approximations to effectively deal with the problems 
associated with spatial and temporal grid refinements. This computer 
program documentation is oriented primarily to someone with this 

level of expertise. In addition to the software description, the 
reference material essential to an understanding of the simulation is 
included in appendices. The first appendix deals with the mathematical 
formulation, the second appendix deals with the numerical methods used 
and the third provides a listing of the computer program with a sample 
run. 


Every computer simulation incorporates simplifying assumptions. 
Complex physical and biological processes must often be lumped into one 
simple empirical equation as for example is done here with root water 
uptake (Hayhoe, 1981). The mathematical complexity involved in 
describing water flow in a layered soil necessitates numerical approxi- 
mations to the mathematical solutions (Hayhoe, 1978a; Hayhoe, 1978b). 
For testing and future development it is important to document the 
details of the simplifications and assumptions that constitute a 
simulation model. The computer software presented here contains a 
unique set of simplifications and assumptions which it is felt merit 
wider application and further development. This technical bulletin 
is viewed as a necessary component in achieving this. 


COMPUTER PROGRAM DESCRIPTION 


The simulation model for soil water movement and uptake by plant 
roots is written in FORTRAN IV using single precision arithmetic. It 
consists of a main program and ten subroutines. The computer program 
is written to be sufficiently versatile that it can simulate a range 
of realistic crop and soil conditions using inputs which one could 
anticipate being available. The input requirements include soil 
description, crop characteristics, weather data and computer program 
controls including spatial and temporal grid refinements. A soil 
water suction head function and a hydraulic conductivity function must 
be specified for each soil layer up to a maximum of four. This can be 
done by selecting functions presented by Clapp and Hornberger (1978) 
based on texture and clay content or by specifying an alternative set 
of functions and parameters which more appropriately represent a 
specific soil. Crop characteristics are input for each crop stage. 
Crop percentage ground cover and root dry matter distribution are 
required. Daily precipitation and potential evapotranspiration must 
be specified. Considerable flexibility has been built in to facilitate 
a range of formats in weather data specification. The output is designed 
to meet the requirements of a soil physicist or an agronomist. A 
detailed output is provided on a daily basis and summaries may be output 
at each iterations as specified by program control. The output 
includes a specification of the total water in the profile, the soil 
water distribution, the plant available water in each layer and the 


estimated suction head, conductivity and diffusivity throughout the 
profile. In addition, the daily and cumulative actual evapotranspira- 
tion, actual transpiration, and infiltration of precipitation are 
output. The daily flux of water through the bottom of the profile is 
output. 


In the following, the main program and the ten subroutines are 
described in sufficient detail to facilitate the set up and running 
of the simulation model as well as its adaptation and further 
development. The equations specified in Appendices A and B are 
referenced throughout by number. 


Main Program 


The main program and the subroutines read the input data on 
logical unit 5. First, the parameters, N, NSTAGE, MD, NLAYR and 
Dit are fread (Format, 415, F5.1) where: 


N is the number of grid-subintervals in the profile; 
NSTAGE is the number of defined crop stages; 

MD is the total number of days to be simulated; 

NLAYR is the number of distinct soils in the profile; and 
DTI is the specified time step. 

Following this ISTAGE(I) is read (Format 515) where ISTAGE(I) is 
the length in days of the Ith crop stage. Four control parameters are 
then read, NOUT, JWAT, IC2, and IC3 (Format 415) where: 

NOUT controls the number of soil water summaries output per day; 


JWAT indicates whether a water table depth is specified - if JWAT 
is greater than zero, water table depths are read; 


IC2 controls an adjustment for water balance - if IC2 is greater 
than zero, no adjustment is made; and 


rc3 indicates the bottom boundary condition - if IC3 equals zero, 
then the content is specified and the flux calculated and if 
IC3 equals 1 then the flux is specified and the content is 
calculated. 


Next the thickness of grid-subintervals in cm, R(1I), is read with 
Format (16F5.2). Following this a card is read for each soil layer 
wreenerormat (12, 6F6.2, 12) where: 


iLe is a parameter used to control the selection or entry of the 
coefficients for the Clapp and Hornberger equations or for 
user specified functions; 


F ; : 3 
FC is the water content at field capacity in cm /cm 

‘ paaeeel: : : 3 3 
PWP is the water content at the permanent wilting point in cm /cm 
WI corresponds to 65 in eq. (4)§ 


DPTH is the thickness of the layer in cm; 


aig 


WK corresponds to 6, in eq. (7) and 


k 


NLAY is the index of the node at the top of the layer described. 


This information is used as input in the subroutine PARM, which is 
now called. 


In order to provide for an examination of the functions resulting 
from the specification in PARM the suction head, conductivity and 
diffusivity are calculated using FUNL and output at selected moisture 
contents for each layer specification. Following this the node depths 
are calculated and output on logical output unit 7. 


Throughout the main program CO7 denotes air dry volumetric soil 
water content. This is set equal to 0.04 cm /cm>. This value could vary 
from soil to soil and in addition one might want to make it a function 
of environmental conditions. The two Gaussian weighting factors, GP, 
are set for future numerical integrations. 


The computer program is written with the objective of being able 
to simulate soil water conditions for a growing crop. For this reason 
the DO loops over time are nested in such a way as to facilitate reading 
the appropriate crop parameters and weather variables. The first loop 
is over crop stages where new crop and root parameters are set. The 
next loop inside the crop stage loop is over days. Here daily weather 
data and daily soil water flux are initialized and evaluated. Within 
the daily loop is a final loop over specified time steps required for 
the numerical solution of the differential equation. Within this loop 
period values for the potential evaporation and precipitation are 
determined in the subroutine PETFIT and the matrix equations specified 
in eq. (15) and (16) are set up and solved for a new soil water distri- 
bution. 


In order to point out some of the special features, more detailed 
aspects of the program will be discussed. Before any timg logps are 
initiated the initial volumetric soil water content in cm /cm~ is read 
(Format 13F5.3) and also printed out on logical output unit 7. In 
addition, the global mass matrix, which is matrix B in eq. (15) is 
set up. Note that since the matrices A and B in eq. (15) are symmetrical 
and tridiagonal only the diagonal and upper diagonal elements are 
calculated and stored. 


At the beginning of each daily increment loop variables are 
initialized as required. Note that in this version of the program 
sunrise in hours, denoted by SUNR, is set equal to 6 and sunset, 
denoted by SUNS, is set equal to 18. For some simulations it would 
be more appropriate to either read or calculate and store actual values 
in the subroutine METDAT for specification here and use in the subroutine 
PETFIT. In addition, the start and termination of daily precipitation 
in hours are denoted with WST and WEND respectively. Here WST is set 
equal to zero and WEND equal to 24. Again if more detailed information 
were available, it could be used. 
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A feature has been added to increase or decrease the size of the 
time step, DT, based on near surface soil conditions and weather 
variables. If less than 0.5 cm of rain occurs and the average water 
content at the top three nodes, CNM, is less than or equal field 
capacity then the time step is doubled. If more than 0.5 cm of rain 
occurs and CNM is less than field capacity or less than the permanent 
wilting point the time step is cut by one half or one quarter respectively. 
This adjusted time step, denoted by DT, is used for the one day. Given 
DT the number of time steps for the day, denote by MP, is determined. 


With each time step eq. (15) is set up and solved and then eq. (16) 
is set up and solved. This is controlled by the index KPC. In this 
version of the program only one corrector step is made but additional 
iterations through eq. (16) are possible by increasing the limit on 
KPC. In addition, it is possible to eliminate the corrector step by 
reducing the limit on KPC required for branching out of the loop. In 
all cases the program sets up the matrix A and S using DP as the 
argument, in eq. (15) or (16). Then eq. (15) or (16) are assembled 
and solved for the new water content at time J+l, CS, using tridiagonal 
Gaussian elimination. The subroutine BCON is then called to estimate 
the water content, CS(1), at the surface. At the bottom boundary for 
the prediction step the water content is set equal to DP(N+1) where 
DP is the water content at time J+3. Alternatively, if the zero 
gradient is assumed, CS(N+l) can be set equal to CS(N). After the 
water content, CS, has been estimated for all N-l nodes, the values are 
checked and forced to lie between saturation and air dry. The index 
KPC is checked. If the corrector steps are complete or not required 
then the time step is completed and KPC is initialized to zero. Other- 
wise, the new estimate for DP is set and KPC is incremented by one. 

Then a new estimate for the water content, CS, is made using eq. (16). 

To complete the calculations associated with each time step CN, the 

water content at new value of J, is set equal to CS and a new DP is 
estimated. The surface actual evaporation, SAE, surface water 
infiltration, SWIN, and the total actual evapotranspiration, including 
root water uptake, TRU, for the period are added to form the corresponding 
daily totals, TSAE, TSWIN and RAE. The control parameter IC3 is checked. 
If IC3 equals 1, CN(N+1) is set equal to the previously specified water 
content, BBC, and the flux through the bottom, FLBP, is estimated based 
Sieene eradient relationship. If IC3 equals 2, then the flux FLB is 

set equal to the specified BBC and the content, CN(N+1), is estimated 
based on the predicted flux. The accumulated flux for the period, FLBP 
and the total for the day TFLB are also calculated. At the end of the 
time step the water in the profile is calculated. An option, whose 

use is controlled by IC2 can be used to check the change in water 

content and force it to match the flux through the upper surface. 
Finally, an accumulation is made of surface flux for the total simulation 
period and a check is made for output of the soil water information for 
the period. When the period iterations for the day are complete then a 
daily summary is output on logical output unit 7. 
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SUBROUTINE PARM 


This subroutine is called once to establish the parameters for 
eq. (4) through eq. (9). Since the program is designed for up to four 
distinct soil layers, there is a loop in the routine which generates 
these coefficients for each distinct layer and stores them in the two- 
dimensional array denoted by CO. The notation for the coefficients for 
each layer are as follows: CO, denotes b in eq. (4), CO, denotes y 
in eq. (4), CO. denotes. 0 “in eq:.).G)>_—C0),. equals ike aim eq. a)i, cos 
denotes c in eq. (6), CO Sdenotes d in eq. (6) and Go_ denotes K. in 
eq. (8)... If a different function ispeci filcatdon ior wre) and K(05 is 
used then up to eight coefficients can be read and stored in CO for 
each distinct layer. In addition, WI denotes 6, in eq. (6) and WKI 
denotes 8, in eq. (7). Calculated values for the coefficients from 
Clapp and Hornberger (1978) are stored in the routine and may be 
selected by specifying the control, ITYPE. Alternatively, if LITYPE 
equals 12, co,» CO, , cO., CO, are read (Format 4F9.6). 


The objective in this routine was to use available information in 
the specification of soil water characteristics. Therefore, if plant 
available water in the Ith layer PAW(1), is specified greater than 
zero, CO, in eq. (4) is adjusted to fit. If the permanent wilting point 
in the ich layer, PWP(I) and field capacity, FC(1), are specified 
preater than zero), then CO. “and CO.cin jeq., @))are cadjustedsto ice 
through these points. Note that the suction at field capacity, SFC, 
and at the permanent wilting point, SPWP, are specified as 15000 and 
344 cm of water respectively. Both SFC and SPWP can be adjusted to 
better fit a particular soil. In addition, where required, the 
coefficient CO, is calculated as specified in eq. (7) and eq. (8). 
une ss the conductivity at field capacity, FCK is specified as 107? cm 
min 


Finally this subroutine calculates as required and prints on 
logical unit 7, the layer thickness, available water, field capacity 
and permanent wilting point for each layer based on the Clapp and 
Hornberger (1978) equations. In addition, the coefficients for eq. (4) 
through eq. (9) for each layer are printed on logical unit 7. 


SUBROUTINE WTABLE 


This subroutine reads the water table depth in cm, DTABLE(I), 
for each day of the simulation. Given the water table depth on day I, 
the program determines the index of the first node lying within the 
water table and stores the result in DTABLE(1). This subroutine is 
called once. If the control parameter, JWAT, is less than or equal 
to zero, no water table is assumed and the subroutine is bypassed. 
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SUBROUTINE METDAT 


This subroutine is called once and its function is to input the 
meteorological data. The present version is structured with the view 
that the basic meteorological parameters could be specified in various 
formats from a number of sources. A control parameter IMET is defined 
and used in a GO TO statement so that routines can be incorporated to 
read data and calculate the required potential evapotranspiration. 

The present version simply reads on logical unit 5, the potential 
evapotranspiration, PEVAP(I), and precipitation, PRECIP(1), for day I 
with format (2F5.2). 


SUBROUTINE ROOTDN 


This subroutine is called at the beginning of each specified crop 
stage and its function is to generate the root density, RDF(I), in 
gm em-3 for the Ith Gauss point. The location of the Ith Gauss point 
is denoted by THG(I), and is calculated in the routine using the grid 
subinterval width, R(1I), and the weighting factor, GP(1). The control 
parameter, IROOT, the fraction ground cover, COVER, and a root effectiveness 
parameter, ROOTEF, are read (Format 13, 2F4.2). IROOT, is used in a GO 
TO statement to select the type of specification of the root. If IROOT 
equals 1, the coefficients are read and the distribution generated with 
an exponential function (Feddes et al., 1974). If IROOT equals 2, the 
root density and depth are read with (Format 2F10.2) and the subroutine 
INTERP is called to generate the root densities at the Gauss points. 
The results are output on logical output unit 7. 


SUBROUTINE PETFIT 


This subroutine is called at each time step and its function is 
to determine the potential surface flux, PET, in cm min ~ for time 
interval Jl. PET is positive for evaporation and negative for precipi- 
tation. The daily potential evaporation, PE, is distributed between 
sunrise, SUNR, and sunset, SUNS, specified in hours using a sine curve. 
The rainfall, W, in cm is distributed uniformly from the start, WST 
hours, and the end, WEND of the precipitation. Daily precipitation is 
distributed uniformly throughout the day. The precipitation for period, 
Jl is then subtracted from the potential evaporation to produce PET. 


SUBROUTINE ROOTUP 


This subroutine is called at each iteration where both PET and the 
fraction ground cover are greater than zero. It calculates the root 
uptake using the model specified in eq. (3) and described in detail in 
the section on mathematical methods. The root suction, ~_, is denoted 
by RP and set to 15000 cm water. Referring to eq. (3), c’is denoted in 
the subroutine by ROOTEF, HKA(I1) denotes the Ith Gauss point value of 
K(®@); HWA(I1) denotes the Gauss point value of (6); RDF(I1) denotes 
the Gauss point value of W(x); and ZAE(11) denotes the Gauss point value 
of R. Note that ZAE is always set greater than or equal to zero and is 
set equal to zero when HWA is less than 50 cm water to account for the 
fact that roots generally do not take up water when oxygen becomes 
limiting. The potential transpiration, X, is determined by multiplying 


PET by COVER. The root uptake, UT, is accumulated until it equals X. 
On the corrector step of the predictor corrector routine, the root 
water uptake in cm for the profile in time DT, UT3, and the cumulative 
root water uptake, TRU, are set. 


SUBROUTINE FUNL 


This subroutine is called once for each iteration. Its function 
is to estimate the hydraulic conductivity, HKA, the soil water suction, 
HWA, and the soil water diffusivity DIA at each Gauss point using the 
Clapp and Hornberger functions parameterized in PARM. In addition, it 
is indicated in the routine where and how one could specify alternate 
functions. 


The subroutine is structured to admit a layered soil. The program 
loops over the number of layers, NLAYR. Firstly it sets the coefficients 
and parameters for the layer and then sets the limits, NMIN and NMAX, 
for a loop over the Gauss points within in the layer. The assumption 
is made that the soil water parameters change is continuous over the 
interface layer. Therefore, the function evaluations at the first Gauss 
point in the new layer are based on the layer above while at the second 
Gauss point on the new layer. This is controlled by the index JMIN. 


SUBROUTINE DISCON 


This subroutine is called once for each iteration if the soil is 
layered. Its function is to add a correction term to the numerical 
approximation to eq. (1) at the soil interface to account for the 
change in soil water suction due to a change in soil. This adjustment 
is required because ) is a function of 9 and x rather than 9 alone 
across the interface layer. The derivative of ~ with respect to x, 
at a fixed time t, across the interface layer must include the 
derivative of ~ with respect to x, denoted by DERIV, for a fixed water 
content, denoted by XAV, as well as the usual derivative of ~ with 
respect to 8 which is incorporated in the diffusivity term in eq. (1). 
The derivative of ~ with respect to x, DERIV, is evaluated across the 
interface layer by assuming linear interpolation between the two layers. 


SUBROUTINE BCON 


This subroutine is called at each iteration. Its function is to use 
the potential flux of water, which is denoted by PET, to calculate both 
the new water concentration at the surface node and the actual soil 
water flux at the surface. If evapotranspiration can occur, the index 
for crop ground cover, COVER, is used to partition evaporation from 
the soil surface and transpiration from the crop canopy. 


Two approaches to estimating the surface boundary conditions have 
been incorporated and they are selected using the control parameters, 
ISURF. If ISURF equals 0 then the equation for water transfer in the 
surface layer is used directly. The concentration at the surface node, 
CS(1), which would be required to give the gradient to produce the 
potential flux is estimated first. Note that HI denotes a weighted 
average hydraulic conductivity and DI denotes a weighted average 
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diffusivity for the surface layer. The weighting is set to best simulate 
evaporation or infiltration conditions. If the resulting value for CS(1) 
does not lie between air dry, CO_, and saturation, CO. then CS(1) is set 
to the appropriate value and flux is estimated based on the resulting 
gradient in water content. 


If ISURF is greater than zero the value of CS(1) is smoothed by 
assuming it is equal to the average water content in the top half of 
the surface layer. The amount of water in the upper half layer is 
budgeted based on flux through the surface and flux of water across 
the mid-point of the layer. 


SUBROUTINE INTERP 


The subroutine is called by ROOTDN if observed root densities 
are used. In addition, it has potential usefulness for interpolating 
between observations of a range of input information. The observed 
dependent and independent variables are denoted within the subroutine 
by Y and X respectively. N denotes the number of observation points 
and M denotes the number of points to be interpolated. Note that the 
routine uses linear interpolation between observed values. The inter- 
polated dependent and independent values are denoted by CINTRP and 
XINTRP respectively. 


RESULTS AND DISCUSSION 


In this section a typical test run of the modelling program is 
presented. It is assumed that to use this computer program, one 
would start with this test case and then modify it to meet special 
requirements. It consists of a 20-day simulation with a layered soil. 
The top 30 cm is assumed to be a sandy clay loam and the remaining 
90 cm is clay. The initial soil water content is set equal to field 
capacity. The potential evapotranspiration is assumed equal to 0.7 cm 
throughout. On day 11 of the simulation, 3.0 cm of precipitation is 
added to the profile. The crop is assumed to go through two 10-day 
development stages during the simulation. The details of setting up 
this computer run are delineated in the following. 


The simulation was run on a VAX mini-computer which is part of 
Agriculture Canada's AGRONET system. The input used for the sample 
run is presented in the appendix. The first statement read from the 
input file sets N equal 18, NSTAGE equal 2, MD equal 20, NLAYR equal 
2 and DT1l equal 60. The next statement specifies the two crop stages, 
ISTAGE(1), equal 10 days each. The third statement sets NOUT equal 3, 
JWAT equal 0, IC2 equal 0 and IC3 equal 0. The next two statements 
read give the 18 grid-subinterval widths denoted by R(1). Two state- 
ments are then read to specify the soil for each layer. For the top 
layer ITP equals 6; FC, PWP and PAW are unspecified and threfore set 
equal 0; WI equals 0.92, DPTH equals 30; WK equals 0.9 and NLAY is 
equal 1. For the second layer ITP equals 11; FC, PWP and PAW are 0; 
WI equals 0.92; DPTH equals 90; WK equals 0.9 and NLAY is equal 8. The 
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next two statements specify the initial values for CN(1). CN is set 
equal to 0.25 in the top 30 cm and 0.373 in the bottom 90 cm which 
correspond to field capacity for the respective layers. Rather than 
reading in the values for precipitation and potential evapotranspiration, 
they are generated in the subroutine METDAT. The final four statements 
are read by the subroutine ROOTDN and specify canopy and rooting characte- 
ristics for the two crop stages. For the first stage two statements 

are read which set IROOT equal 1, ROOTEF equal 1, COVER equal 0.5, A 
equal 5.13 and B equal 0.32. The final two statements in the input file 
specify for the second crop stage, IROOT equal 1, ROOTEF equal 1, COVER 
equal 0.9, A equal 5.13 and B equal 0.08. The coefficients A and B for 
the two crop stages, give rise to root density distributions as plotted 
in figure 1. 


The output generated by the soil water model run with these 
controls is included in the appendix. The first two output lines are 
the input statements specifying the soil and layer thickness. The 
layer characteristics as calculated in the subroutine PARM are then 
printed, including the coefficients for the soil water functions. In 
this simulation field capacity is set to 344 cm of water and the 
permanent wilting point to 15000 cm of water. In addition K, in eq. (9) 
is set equal to 10-2 cm/min. In order to specify more elearty the 
resulting soil water functions for each soil layer, their value at ten 
equally spaced water contents is calculated and output and plotted in 
figures 5, 6 and 7. The node depths and initial water content at this 
node are printed. 


The root distribution generated in the subroutine ROOTDN assuming 
an exponential function is printed for the first crop stage. This is 
followed by detailed results of fluxes, root water uptake and soil 
water distribution estimation. Since NOUT is set equal to 3, the 
estimated water content at the nodes is output as well as the water in 
the profile and the flux through the bottom every eight hours. Note 
that in this sample model run, zero water content gradient is assumed 
across the lowest grid subinterval. It should be pointed out that since 
the computer simulation involves a numerical approximation to eq. (1), 
that inappropriately large spatial or temporal grid refinements can 
lead to unacceptable errors. In this output hourly time steps were 
used and the errors resulting from the numerical techniques did not 
appear unacceptable given the levels of accuracy expected in the input 
data. For large hydraulic conductivities, used with large time steps, 
the estimated flux through the bottom can be too large and should be 
checked against the change in the water in the profile. 


At the end of each day a summary of the fluxes and soil water 
profile is printed. Firstly the estimated average water content and 
the percent plant available water are printed for each soil layer. 
A summary of daily and cumulative values for potential evapotranspira- 
tion (PE) and actual evapotranspiration (AE), transpiration, precipitation 
and infiltration is printed. 
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Since crop stage 2 starts on day 1l, the new crop root distribution 
is evaluated and printed. For stage 2, the crop canopy parameter, 
COVER, is set equal to 0.9. In the flux summary for day 11, the results 
of the first precipitation are printed. The precipitation is 3 cm and 
2.3 cm infiltrated the soil. The difference equals the PE which is 
subtracted before the water enters the soil. Under these conditions 
the model leaves the AE and transpiration remain at 0 for the day. 


The accumulated volumes of potential evapotranspiration, actual 
evapotranspiration and its components are shown in figure 2. Although 
during the first 2 or 3 days of the simulation run the actual evapora- 
tion and the actual transpiration were very similar, later on the rate 
of actual evaporation dropped off quickly to less than 0.1 cm day 
Even the 3 cm rainfall on day 11 did not substantially increase this 
rate. Consequently, the rate of actual transpiration approximated 
the rate of actual evapotranspiration. 


During the first 10 days of the simulation period, the rate of 
actual evapotranspiration decreased slowly (Figure 2) in response to 
decreasing water availability (Figure 3), especially in the upper layer. 
However, because of the increased rooting depth and density and the 
3 cm of rain on day 11, the rate of actual evapotranspiration followed 
closely the potential evapotranspiration rate. For example, the model 
predicted an actual average evapotranspiration rate of 0.65 cm day 1 
between day 11 and 20, while the potential rate was 0.7 cm day. 


Some of the soil water content profiles produced by this particular 
simulation run are shown in figure 4a and b. The discontinuous nature 
of these profiles is due to the fact that we assumed that the soil water 
potential and not water content distribution itself, was continuous 
throughout the profile, across layer interfaces. During the first 10 
days the water content in the upper layer decreased rapidly by 72% 
to an average water content of 0.18 em? cm73, while that of the lower 
layer decreased by only 2% to 0.35 cm3 cm73. Plant available water 
dropped correspondingly from 100% to 33% in the upper horizon and to 
73% in the lower layer (Fig. 3). During the second half of the 
simulation run, when the root system had expanded SESS IG a the water 
content of the upper layer decreased to 0.13 cm? cm » which was below 
the permanent wilting for the layer, while the water content of the 
lower layer had decreased to 0.32 cm? cm ~, which corresponded to 54% 
plant available water. The water in the entire profile fell from 37.6 cm 
on day 10 to 34.1 cm on day 20. However, the ratio of actual to 
potential transpiration increased from 0.65 to 0.91 during'i the same time 
period, which demonstrates the model's sensitivity to changes in the 
root system. 


DEPTH IN cm 


80 


100 


Figure 1: 


ROOT DENSITY (g cm-3) 


@ STAGE 1 
@ STAGE 2 


Root density distributions used in the simulation to 
demonstrated root water uptake 


CUMULATIVE VALUES (cm) 


14 


2 


10 


Figure 2: Simulated evaporation and transpiration for 20-day period 
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Figure 3: Simulated plant available water over a 20-day period for both the 
upper and lower soil horizon. 
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Figure 4a: Simulated soil water content profiles at selected 
times during the first 10-days 
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Figure 4b: Simulated soil water content profiles at selected 
times during the last 10-days 
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Figure 5: Soil water and suction head relations 
used in simulation. 
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Figure 6: Hydraulic conductivity functions used in 


simulation 
SANDY 
GRAN, 
LOAM GIPAN, 
0.1 (072 0.3 0.4 0.5 


VOLUMETRIC WATER CONTENT (cm3cm-3) 


0.6 


DIFFUSIVITY ( cm2min-") 


10 


109 


10-1 


10-2 


10-3 


1034 


10-5 


Figure 7: Diffusivity functions used in the 
simulation 
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APPENDIX A 
Mathematical Methods 


In this section we provide a summary of the mathematical procedures 
used in the modelling program. The mathematical equations for the 
movement of soil water and the uptake by plant root systems have been 
specified as follows: 


28 8 (D(6) 36 + K(®)) + R(@, x, t) (1) 
dt we ROX ox 
where 
t is time in min, t takes values from 0 to T; 
6 is the volumetric water content in cm? of water per ae of ssoils 
x is the vertical distance from the soil surface in cm (positive 
downward), x takes values from 0 to L; 
y(8) is the soil-water suction head function in cm; 
K(6) is the hydraulic conductivity of the soil in cm per min; 
D(8) is the soil-water diffusivity function defined by the equation 


dy ‘ 

D(6) = Go K(8) (2) 
R(6, x, t) is the water uptake by roots in ae of water per me of 
soil per min. To solve eq. (1), 9(x,t) must be specified for all x at 
t=0 and for x=0 and x=L for t=0 to T. Meteorological conditions are 
required to specify 8(0, t) and L is chosen such that an appropriate 
value for 68(L, t) can be determined, either by specifying a concentration 
or flux condition at x=L. Following Feddes et al. (1974), the root sink 
term is defined by: 


R(9, x, t) = cK(®) (W - ¥(6)) WO) (3) 


where 3 
W(x) is the root weight function in gm per cm ; 
yw. is the root suction head; 

c is a parameter to adjust for root effectiveness. 


A major obstacle to successful application of this mathematical model 
to an actual soil system relates to the problem of appropriately specifying 
K(8), ~(8), and R. A potentially useful approach to this problem has been 
suggested by Clapp and Hornberger (1978). They propose standardizing 
K(6), and ~(6) so that a change in soil would only require a change of the 
parameters in the standardized relations. Further, they provide parameters 
for their functions for a wide range of soil textures and soil clay 
contents. The problem of specifying K(®@) and (8) now reduces to that of 
providing the appropriate parameters by selecting them from a table 
according to textural group (Clapp and Hornberger, 1978) or by deriving 
them from point observations on K(6) and (8). Following Clapp and 
Hornberger (1978), 


-b 
Wy (C) = Ee (eye) O= 6 <6 (4) 


2b+3 


K(6) = K, (0/8,) (5) 
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where b,. 0 ,.¥ .~and Ke are the parameters to be specified. It is assumed 
that 6 is the saturated water content and that K_ is the saturated 
hydraulic conductivity. Since ~(§_) = 0, a value’6., is selected near 

6 _, and a parabola for (6) is fit°for @ between 9, and 6 such that (6) 
and its derivative are continuous at O.. On this interval, 


y(@) = c (6/8. - d) (8/6. - 1) EG SoGe< 8G (6) 
where c and d are evaluated according to the conditions just specified. 


During the course of this study, it was found that a similar 
adjustment might be required for K(6). For example, if K(6) undergoes 
a very sharp increase near 6_, as is commonly observed in highly 
structured media and fine textured soils, selecting K_ equal to the 
observed saturated hydraulic conductivity would resulf in values for 
K(6) which are too large for a wide range of intermediate 6 values. 

The adjusted equation for hydraulic conductivity was then specified by: 


K(6) =°K, (8/0,)°>*9 Ons 0 <6, (7) 
K(8) = Ky (0/0) "> + (K, - Ke) ((6-6,)/(8, - 8,))° 
ae Sas ae 


where 6, is an intermediate value for 9 for which K(6, ) is known or 
can be estimated. The value of the parameter Kes is doce emned from 
eq. (7) by the relation 


2b+3 
Ke = K(6,) Cy Cy) (9) 


For this study, 8, was selected such that (6,) equalled 344 cm 
H,0, a value which has been found to approximate field capacity (Colman, 
1947; Haise et al., 1955). Field capacity was chosen because by defini- 
tion it corresponds to a condition of low conductivity for all soils. 
Data presented by Ahuja (1973) and Bruce_and Whisler (1973) suggest 
that this small value may be between 107° and 10-© cm min7! for a wide 
range of soil textures. Hence, in the absence of a measured value for 
K(8,) it was assumed to be ees cm min ~, except under the condition 
where the resulting K, was greater than Kos in that case eq. (5) was used. 


It should be pointed out that the estimated values for (6,) and 
K(6, ) are averages. Undoubtedly better estimates can be obtained and 
substituted when specific soil types and textural groups are considered, 
and of course measured data should be used whenever they are available. 


A major function of the mathematical model for soil water is to 
predict the actual evapotranspiration. The potential evapotranspiration 
is the combined evaporation from the soil surface and transpiration from 
the plant which would occur if water were not limiting. It is determined 
by meteorological and crop parameters (Feddes et al., 1974). The 
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partitioning of the potential evaporation between the plant and soil 
surface is determined by the crop cover similar to the procedure 
proposed by Nimah and Hanks (1973). 


Water is assumed to evaporate from the soil surface at the 
potential rate when moisture is not limiting which provides a flux 
surface boundary condition. Under drier conditions the actual evaporation 
is controlled by the rate of upward water movement to the soil surface. 
Under such conditions, the upper boundary condition is to set the soil 
water content to its air dry value. Similarly, for precipitation water 
enters the soil at a rate equal to the precipitation until the surface 
becomes saturated. Then the surface is assumed saturated and infiltra- 
tion is limited by the rate with which water flows through the surface 
zones. 


Transpiration by the plant is controlled by the sink term specified 
by eq. (3). Feddes et al. (1974) proposed using an iterative technique 
to evaluate y_ such that the plant would transpire at the potential rate 
until y reached a maximum of 15000 cm of water. The assumptions 
specified in their model appear to permit water moving from the plant 
into the soil. Here, a simpler procedure is proposed. It is assumed 
that the crop will attempt to meet evaporative demand from the top zone. 
If this is not possible it will extract water successively through lower 
zones until the potential transpiration is achieved. If the potential 
is not achieved then the actual transpiration is the sum of the water 
extracted from each zone. The model assumes the value of ~ in eq. (3) 
reaches the maximum of 15000 cm of water in each zone. This model speci- 
fication should not be interpreted to mean that all the plant available 
water is extracted from the upper zone before any water is extracted 
from lower zones. Firstly, whenever the actual transpiration is less 
than the potential, water is extracted from all zones simultaneously as 
specified by eq. 3. These results would be identical to those predicted 
by the model presented by Feddes et al. (1974). When water is available 
both near the surface and deep in the profile, the model assumptions 
imply that the crop would first use the readily available moisture near 
the surface. The justification for this assumption is that in the top 
layer of the soil the temperature, nutrients and aeration are generally 
most favorable for root activity. As water is depleted, uptake from the 
surface zones decreases rapidly due to the rapid decrease in K(8) in 
eq. (3). The result is again a more uniform extraction pattern of soil 
water by the crop roots. 
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APPENDIX B 


Numerical Methods 


Eq. (1) is solved numerically using linear piecewise elements and 
the Galerkin technique as described by Hayhoe (1978). The space variable 
is discretized by using linear elements and the time variable is 
discretized using a Crank-Nicolson procedure. 


Consider first the general approach. Let S be the set of real- 
valued, piecewise continuously differentiable functions, w(x), on 
(0,2), and let S_ denote an m-dimensional subspace of S which is 
spanned by m basis functions, w.(x). Assume a solution to eq. (1) of 
the form cs 


m 
oP (x,t) = a a, (t) wy Ox) (10) 


The a, (t)'s are determined in part by the initial and boundary conditions. 
An operator, L, is defined on 6 by 


L(o) = 22-2 (p(o) 32 + K(e)) - R(O, x, t) (11) 


The remaining a. (t)'s are determined from the relationship 


(L(6), i) = J L(6 ) w, dx = Okita = 2 Smeets tin i C2) 


for all o<t<T. Rearranging eq. (12) and integrating by parts one obtains 
a system of m-2 first-order differential equations which in matrix 
notation can be written as: 


da(t) 


P dt 


+ A(6_) a(t) = S(6_). (13) 


Here B is an (m2) x (m-2) matrix with (i,j)th element equal 


y 
i “Ses A(8_) is an (m-2) x (m-2) matrix with 
71 
(i,j)th element equal to 9 dw. dw. 
ji DiCe: ) —1 dx; 
m dx dx 
fe) 
a, (t) 
a, (t) 
a(t) is an (m-2) x 1 matrix equal to . sand S(6_) is an 
a (t 
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(m-2) x 1 matrix with jth element 


I da, Ct) rear CE) “a oh + D(6) Ge en? oe 
A dt dt dx dx 


OY eox( eee tio=In (eee ke We ae 
dx 5 dx , J 


For ‘thesmatrices.. inveqwrG3)2is=7293) a22,m—! Vand j= 243,.50,0- i 


For linear elements, chapeau functions are used. The interval 
(o,2) is divided into N subintervals with h, denoting the length of 
the ith interval. The functions, w., are chosen so as to form what 
is called a local or patch basis since they are zero everywhere in 
(o,2) except on at most two adjacent subintervals. The basis function 
is defined by the relation 


i-1l 
Wiyd (x) = (xa=a9 1S hy) i h, on the ith subinterval (14) 
k=1 
i+l 
(tee het=ex) babe on the (i-1)th subinterval 
kel i i+] 


0 elsewhere on (0,2). 


Substitute these basis functions into eq. (13), A(O_) and B become 
(m-2) x (m-2) symmetric tridiagonal matrices. 


The next step in determining the numerical approximation to eq. (1) 
involves the discretization of eq. (13) with respect to time. This was 
achieved using Crank-Nicolson extrapolation (Douglas and Dupont, 1970). 


Define 
j LJ : = 1S Qi >m 
a, = a, Gat) EO ae 
j 
and = ‘ 
n | (x) = ce (x, jAt) 
where At = T/K. Then (13) becomes, 
; f : x alia {-2 j 
e a ob ELS GS oS a) ate =| 5 Be) A a 
a m m = At 
At 
alle ote | 
+ S(1.5 on A) oe ) (1: 


An optional corrector step can be indicated using the estimate from 
eq. (15), which is denoted by gJ* (Douglas and Dupont, 1970). The 
corrector step is specified by 


RG ies 


(F Hae 5 A(.5 (el + 3*))) alti 
(16) 

‘ ne ; , ei 

“(52 B= 5 A(.5 (G2 + git) a + 5(.5 (@) + ery) 

The problem of using the finite element method to find a numerical 

solution to eq. (1) has now been reduced to solving a linear matrix 

equation for each j. When j = 0, a is determined using nodal values 

from the initial condition and thus 6° is known from eq. (10). Since 

ed and ei-4 are required in eq. (15) ™ to determine ees, eq. (G5) 

cannot be applied to the calculation al, For this initial case a 

predictor-corrector technique is used (Douglas and Dupont, 1970). For 

aeeis2s--~> KL, eq. (15) and eq. (16) define the procedure for 

evaluating as » and from the boundary conditions al and ad are evaluated. 


An obstacle to the efficient application of Galerkin procedures 
for nonlinear equations of the type presented here is that the matrices 
in eq. (15) must be calculated at each time step. This involves 
performing an integration such as 


x ’ : 
a Dee cee see) Gx (17) 
; m m == === 
fe) dx dx 
for all i and k equal 2,3,... m-1. We considered this problem and 


tried to make our solution as efficient as possible. The integrals 
were calculated using Gaussian integration (Dahlquist et al., 1974) 
with two Gauss points per subinterval for linear elements. The result 
of this was that for example, only D(1.5 Euper <5 als ) had to be 
calculated for each Gauss point and then the required summations 
performed to evaluate A(8_) at each time step. 


As has been pointed out, most of the elements in A(O_) and B are 
zero. This was taken advantage of in order to economize on program 
size and to improve efficiency. A technique described by Schwarz et 
al. (1973) for banded matrices was used to solve eq. (15) and eq. (16). 
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